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ABSTRACT 



Context. Very high-energy (VHE, E > 100 GeV) y-ray data are a valuable input for multi-wavelength and multi-messenger (e.g. com- 
bination with neutrino data) studies. 

Aims. We aim at the conservation and homogenization of historical, current, and future VHE y-ray-data on active galactic nuclei 
(AGN). 

Methods. We have collected lightcurve data taken by major VHE experiments since 1991 and combined them into long-term 
lightcurves for several AGN, and now provide our collected datasets for further use. Due to the lack of common data formats in 
VHE y-ray astronomy, we have defined relevant datafields to be stored in standard data formats. The time variability of the combined 
VHE lightcurve data was investigated, and correlation with archival X-ray data collected by RXTE/ ASM tested. 
Results. The combination of data on the prominent blazar Mrk421 from different experiments yields a lightcurve spanning more 
than a decade. From this combined dataset we derive an integral baseline flux from Mrk421 that must be lower than 33 % of the Crab 
Nebula flux above 1 TeV. The analysis of the time variability yields log-normal flux variations in the VHE-data on Mrk421. 
Conclusions. Existing VHE data contain valuable information concerning the variability of AGN and can be an important in- 
gredient for multi-wavelength or multi-messenger studies. In the future, upcoming and planned experiments will provide more 
data from many transient objects, and the interaction of VHE astronomy with classical astronomy will intensify. In this con- 
text a unified and exchangeable data format will become increasingly important. Our data collection is available at the url: 
http : //nuastro-zeuthen. desy . de/magic_experiment/pro jects/light_curve_archive/index_eng . html. 

Key words. Gamma-rays: observations - Galaxies: active - Galaxies: individual: Mrk421 - Galaxies: individual: Mrk501 



O 1. Introduction 



The broad-band emission observed from active galactic nuclei 
(AGN) spans the complete electromagnetic spectrum from radio 
to VHE (very high-energy, E> 100 GeV) y-rays. Since the dis- 
covery of the blazar Mrk421 in the VHE regime (Punch et al. 
1992), many new detections of AGNs by different experiments 
have been reported (see, e.g. Horns 2008; Hinton & Hofmann 
2009, for reviews). Strong flux variations on different timescales 
were observed from many AGNs. Variability in the VHE y-ray 
regime was measured down to the minute timescale (MAGIC 
Collaboration et al. 2008; Aharonian et al. 2007). Many AGN 
lightcurves were produced by different y-ray experiments from 
observations of flaring states, dedicated monitoring programs, or 
from joint campaigns with other experiments. Long-term y-ray 
lightcurves are valuable to address open questions concerning 
AGNs, such as possible periodicities (see Thieler et al. 2010, an 
analysis using the data collection of the present publication), the 
log-normality of their flux state distribution, the nature of the 
radiation mechanism of AGNs, or the estimation of their duty- 
cycle. 

A log-normal behavior can be indicative of a multiplica- 
tive process of the underlying mechanism governing the vari- 
ability of the object. In the case of AGNs, this could be evi- 



dence for a connection of the observed emission to accretion 
disk activity (see Giebels & Degrange 2009, and references 
therein). Previously, log-normal flux variations were reported 
in the X-ray band from two objects: the narrow-line Seyfert 1 
galaxy IRAS 13244-3809 (Gaskell 2004) and BLLac (Giebels 
& Degrange 2009). In VHE data, log-normal variability was ob- 
served from a high flux state of the BLLac object PKS 2155-304 
(Degrange et al. 2008). 

The nature of the AGN radiation mechanism in the VHE 
regime, i.e. whether the observed radiation has a leptonic or 
hadronic origin, still remains ambiguous. The detection of neu- 
trinos from these objects would prove the existence of a hadronic 
component. In the framework of multi-messenger strategies 
(e.g., combining electromagnetic with neutrino data), the phe- 
nomenology of lightcurves in the electromagnetic wavelength 
band can give valuable input. The Neutrino triggered Target 
of Opportunity (NToO) program (Bernardini 2005; Ackermann 
et al. 2008) is based on the idea of neutrino events from a vari- 
able object (single events or multiplets) that are used to trigger 
VHE monitoring of the same object. Coincidences between neu- 
trino triggers and y-ray flux high states, which occur more often 
than expected from random coincidence with atmospheric neu- 
trinos, would be evidence of a hadronic component of the VHE 
y-ray signal and a cosmic origin of the neutrinos. Long-term 
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y-ray lightcurves can be used to estimate the AGN dutycycle, 
which is an important input parameter for such analyses. 



2. Lightcurve data 

Considering the heterogeneous nature of historical and present 
VHE data comprising different file and content formats from 
different experiments, a common data format is desirable. 
Moreover, as described above, VHE astronomy has already 
started to interface strongly with different fields of classical as- 
tronomy. In view of an effective exchange and diffusion of VHE 
data within the astronomy and astroparticle community, our stor- 
age strategy is to use an ascii file format (SLF) and the widely 
used FITS and VOTable file formats (see following section). 

The use of a standard data format will also become in- 
creasingly important in the framework of the next-generation 
Cherenkov telescope systems, such as the upcoming CTA 
(Cherenkov Telescope Array) experiment (see e.g. Martinez 
2008). In the CTA era, VHE astronomy will intensify multi- 
wavelength (and multimessenger) interactions with other fields. 
Furthermore, a standard data format is essential for running an 
experiment like an observatory, and make it open to the whole 
astronomy community. 

Publicly available lightcurve data from 1992 until today were 
collected from the Whipple (Kerrick et al. 1995; Schubnell et al. 
1996; Buckley et al. 1996; Maraschi et al. 1999a), HEGRA 
(Aharonian et al. 1999b,a; Krawczynski et al. 2001b; Aharonian 
et al. 2001, 2002, 2003, 2004; Kestel 2002), CAT (Piron 2000; 
Piron et al. 2001), H.E.S.S. (Aharonian et al. 2005, 2006), 
MAGIC (Albert et al. 2008; Donnarumma et al. 2009), and 
VERITAS (Rebillot et al. 2006; Donnarumma et al. 2009) exper- 
iments. We are also working on collecting more data from these 
and other experiments such as the Crimean GT-48 Cherenkov 
telescope (e.g. Neshpor et al. 2007), the Tibet air-shower ar- 
ray (e.g. Amenomori et al. 2003), or the Patchmarhi Cherenkov 
Telescope Array (e.g. Gupta et al. 2008). 

2.1. Data collection and formatting 

Different formats and standards are used by different experi- 
ments. We use simple directly usable ascii tables and the stan- 
dard astronomy file formats FITS and VOTable to store the col- 
lected, combined data. In the present work, the types and units 
of the FITS/VOTable datafields (i.e. columns) follow a simple 
lightcurve format (hereafter SLF) given in Table 1 . The data are 
also provided in ascii file form following the SLF definition. 
These files are referred to as slf-files. The tables contain rows 
of (typically nightwise) integral flux-measurements and further 
information necessary for subsequent analyses. The units of the 
datafields are defined within the FITS/VOTable files and follow 
the conventions defined in Table 1 within the slf-files (ascii). 



The SLF datafield format definition provides the means to effec- 
tively combine heterogeneous datasets. In case further informa- 
tion than included in columns 1 to 12 are necessary, additional 
entries are foreseen. For example, to define a spectral shape dif- 
ferent from the standard pure powerlaw, the spectral parameter- 
ization together with the additional parameters must be defined 
as comments within the data files. In case of lightcurve data only 
containing information on a few datafields, the flexibility of the 
binary FITS and VOTable file format allows minimizing the data 
volume. For example, if only the observation date and the flux 



value are known for the total lightcurve data of an object, the 
binary file can only contain MJD star , and F, along with the def- 
inition of units used. While this example represents the simplest 
case of a homogeneous dataset, in practice, the heterogeneous 
lightcurve data combinations from different experiments can in- 
clude information on a specific datafield for one subset of the 
combined data (i.e. for one experiment), and not for another sub- 
set. In this (most common) case, when information is not avail- 
able for datafields of some of the measurements, the following 
conventions for default datafield values apply. 

- If no start and end times are known, then MJD star , is set to 
the same value as MJD m id, i.e. the middle of the exposure 
time, and MJD enc i to -1. In this case the duration might still 
be given. 

- In some cases only the MJD of the observations day is 
known. Then, MJD m id is the MJD of the observation day. 

- Most energy spectra can be described by a pure powerlaw 
(differential spectral index a) or a powerlaw with an expo- 
nential cutoff (cutoff energy E cuto ff). Additional fields are 
foreseen for cases with different parameterizations of the en- 
ergy spectrum. If such a different parameterization is given, 
a defining formula containing the relevant fields has to be 
written in the comment section of the lightcurve file. 

- The flux-flag (fflag) is a one-character flag indicating 
whether the integral flux F is an upper limit (fflag = '<') 
or a flux measurement (fflag = '='). 

- In the present work, we chose Crab Nebula flux units for the 
observed integral flux F above the energy threshold E t h r . 

- For any other entry, the value -1 means that no information 
is available. 



2.2. Lightcurve combination 

For a combination of lightcurves from different experiments, the 
measured integral flux values must be converted to a common 
energy threshold. For the results in this paper, we chose an en- 
ergy threshold of 1 TeV The conversion of integral fluxes in Crab 
units to the same energy threshold requires knowledge of the en- 
ergy spectrum of the considered object and of the Crab Nebula as 
measured by the same experiment. The differential energy spec- 
trum of the Crab Nebula in the energy range covered by most 
VHE y-ray experiments is described well by a pure powerlaw 
of the form <f>cmb(E) = (p^lab ' (E /TeVy rcmb in the energy range 
from 100 GeV to few tens of TeV. Deviations from the pure pow- 
erlaw form at both ends of this energy range (as expected in the 
framework of the inverse Compton scenario) were not taken into 
account here. Extrapolation to energies below 100 GeV requires 
a full parameterization of the spectral energy distribution of the 
Crab Nebula (e.g. Aharonian et al. 2004). An integral flux in 
Crab units above the energy threshold E,h r is given by 



F(E > E thr ) = 



C dEc/>(E) C dE<p(E) 



C dE(f>Crab(E) FcrabiE > E thr ) 



(1) 



where 0(E) and <pcmb(E) are the differential energy spec- 
tra of the source and of the Crab Nebula (in units of pho- 
tons TeV _1 cm _2 s _1 ). The corresponding integral fluxes above a 
given energy threshold are denoted with F(E > E,h r ). In the 
case of a pure powerlaw differential energy spectrum (f>(E), the 
conversion of an observed integral flux F(E > Eq) to a given 
threshold energy E,h r is given by 

r+1 F(E > Ea~) 

(2) 



Fit >£,,„.) ={^- 



F C rab(E > Ethr) 
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Dataneld 


Name 


TTVDF 

1 1 irh, 


TT TIVTT'T' 


Explanation 


i 
1 


MJD mid 


mjd_mid_exp 


MJD 


center of observation exposure 


L 


MJD start 


mjd_start 


MJD 


start of observation 


3 


MJD md 


mjd.end 


MJD 


end of observation 


4 


r 


int_flux 


Crab units 


Observed integral flux 


c 

J 


A L~ 

&F stat 


sigma_int_flux_stat 


Crab units 


Statistical error on F 


s- 




A TP 


sigma_int_flux_sys 


Crab units 


Systematical error on F 


1 


a 


alpha 


none 


Differential spectral index 


8 




sigma_alpha_stat 


none 


Statistical error on a 


9 


Aa sys 


sigma_alpha_sys 


none 


Systematical error on a 


10 




e_thr 


TeV 


Energy threshold in TeV 




LP 

^cutoff 


e_cut 


lev 


Cutoff energy in TeV 


12 


Experiment 


experiment 


none 


Experiment string 






UUIaLlUIl 


hours 


UUIaLlUIl Ul UUSClvallUll 


14 


Reference 


reference 


none 


Reference string (pref. ADS format) 


15 


Flux-flag 


mag 


none 


flux-measurement: ' = ', upper limit: '<' 


16 


additional entry 1 


add_entryl 




e.g., special parametrizations 


17 


additional entry 2 


add_entry2 







Table 1. Datafield definition* used for storage of lightcurve data in the FITS/VOTable file format. 

* The used datafield definition is referred to as simple lightcurve format (SLF) in the text. TTYPE and TUNIT are the name and units used within 
the FITS file. 



where Eq is the energy threshold of the observation. 

The combined day-wise, integral flux lightcurve of the 
BLLac object Mrk421 above 1 TeV is shown in Figure 1. This 
lightcurve includes all data we have collected so far, cover- 
ing an uprecedented 17-year time-span from 1992 to 2008. In 
Figure 2, the day-wise integral flux lightcurve of the BL Lac ob- 
ject Mrk501 is shown. 

2.3. Systematic errors 

The inhomogeneity of data combined from different experiments 
and the inclusion of data from the pioneering time of VHE y-ray 
astronomy induces systematic errors that are not easily evalu- 
ated. Combination of data from many different experiments, dif- 
ferent time periods, partly inaccurate flux normalizations (i.e. 
y-ray rate measurements) and energy thresholds can lead to sys- 
tematic errors in the relative flux normalizations. When com- 
bining data, the flux normalization depends on assumptions on 
the spectral shape and the energy threshold. Here, limiting fac- 
tors can be experimental uncertainties (especially in old data) 
induced by poor knowledge of atmospheric conditions, aging de- 
tector effects, and intrinsic spectral variability of the considered 
objects. We convert flux and rate measurements to Crab Nebula 
flux units using Crab Nebula data taken by the same experiment 
as close in time as possible, therewith reducing systematic un- 
certainties induced by seasonal or instrumental variations. 

The error on the reconstructed spectral shape results in an 
error on the flux state when converting to different thresholds. 
For example, when converting integral fluxes from 1 TeV to 
100 GeV, typical systematic errors on the spectral index of 0.1 
(pure powerlaw case) lead to a relative systematic error on the 
integral flux of 30 % due to this extrapolation. The datafields 
contained in our FITS files provide all information needed to es- 
timate these systematic errors for conversion to different energy 
thresholds. A conversion over a narrower energy range than used 
in this example limits this systematic error to less than 30 %. 

The flux units used in some publications are given in un- 
calibrated units (e.g. counts per minute), i.e. as an instrument- 
specific y-ray count-rate. The energy threshold is also unknown 
in some cases. For flux calibration of lightcurves from the 
Whipple experiment before 1998, we used Crab Nebula data 



taken by the same experiment as close in time as possible (see 
Hillas et al. 1998, and references therein). However, due to de- 
pendencies on zenith angle and weather conditions, normaliza- 
tion of rate measurements is very uncertain. We estimate that 
the systematic error induced when including rate measurements 
is below 40 %. This form of relative systematic error can be 
avoided when restricting analysis to flux measurements, by re- 
stricting to data from a single experiment, or by intercalibrating 
the different measurements when observations from different ex- 
periments overlap. 

3. Discussion 

Significant evidence of a correlated gamma- and X-ray emission 
of blazars was presented in earlier studies (see e.g. Takahashi 
et al. 1996; Buckley et al. 1996; Maraschi et al. 1999b; 
Krawczynski et al. 2001a; Neshpor et al. 2007). Usually the re- 
ported correlation is linear, but in a few cases a quadratic re- 
lation between the fluxes in both bands was found (see e.g. 
Krawczynski et al. 2000). Such a correlation provides us with 
essential information on the underlying acceleration and emis- 
sion processes and is especially valuable for variability stud- 
ies. Very often the gamma- X-ray correlation is interpreted as a 
strong argument in favor of the so-called synchrotron-Compton 
jet emission models in which the same population of ultrarel- 
ativistic electrons is responsible for production of both X-rays 
via synchrotron radiation and TeV y-rays via inverse Compton 
scattering (Katarzynski et al. 2005). However, it can also be ac- 
commodated in the hadronic framework in particular in mod- 
els which assume that the observed y-ray emission is a result 
of interactions of accelerated protons and ambient gas or low- 
frequency radiation (Aharonian 2000; Miicke et al. 2003). 

Using contemporanous X-ray data extracted from the 
RXTE/ ASM database web interface at MIT 1 , we calculated cor- 
relation coefficients between the VHE and X-ray bands. For each 
VHE measurement, an average ASM count rate was calculated 
centered on the VHE observation date (± 6 hours). Only using 
measurements with a significance of at least 3 standard devia- 
tions for the VHE data and 10 counts for the ASM measure- 

1 http://xte. mit. edu/ASMJc. html 
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Fig. 1. Long-term lightcurve of Mrk 421 (day- wise integral flux). Data from the major y-ray telescopes were combined and normalized to the 
same energy threshold (1 TeV) and converted to Crab units (see text). A zoom into the period of strong activity (2000/2001) is also shown. 
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Fig. 2. Long-term lightcurve of Mrk 501. Available data were combined and normalized. Shown are the day- wise integral fluxes above 1 TeV in 
units of the Crab Nebula flux. 



ment, these data yield correlation coefficients (Spearman rank) 
of 0.65 for Mrk 421 and 0.68 for Mrk 501. In both cases, the 
probability of an uncorrelated system producing datasets with 
similar Spearman rank correlation coefficients is less than 10~ 8 . 
The flux measurements in both wavelength bands are not strictly 
simultaneous. However, they represent a measure for the aver- 
age daily flux state of the objects, and the observed behavior is 
consistent with correlated average VHE/X-ray daily flux levels 
modulated by shorter term flux variations. 

In Figure 3, the distribution of the VHE flux states of 
Mrk 421 (also see Figure 1) is shown. Each entry of this his- 



togram represents a snapshot flux measurement of Mrk 421 as 
given by the combined dataset. All individual observed inte- 
gral flux values were converted to flux values in units of the 
Crab Nebula flux and normalized to a common energy threshold 
of 1 TeV. The observed energy spectra of the Crab Nebula and 
Mrk 421 as observed by the individual experiments were taken 
into account. The overall distribution can be described by an ex- 
ponential law (left panel) above a flux of a few tenths of Crab. An 
exponential shape might indicate that the measured flux states 
mainly reflect a stochastic outburst state of the object. As shown 
in the right panel of Figure 3, a better fit to the data is obtained 
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Slope 
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°Gauss 
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VHE 



>1TeV [Crab] 



Fig. 3. Distribution of VHE flux states of Mrk421. An exponential function fit above a flux of 0.25 Crab, (avoiding detector threshold effects) can 
describe the data (left). The data are very well described by a fit of a Gaussian+log-normal distribution (right). 



when using the sum of a Gaussian and a log-normal distribution 
(Aitchison & Brown 1957; Limpert et al. 2001), as given by 



tt , N Ln / (log(x)-tf \ 

fix) = =exp — . 

xo-^i \ 2cr 2 / 



(3) 



The data were further divided into time intervals of equal length, 
each comprising 40 flux observations with statistical errors 

cr ; . The excess variance cr xs = ^jjj Tli = \(fi ~ f) 2 _ °f (Vaughan 
et al. 2003) was calculated for each time interval. The ex- 
cess variance cr xs is a measure of the Poisson noise corrected 
rms value of the fluxes within the corresponding time inter- 
val. Figure 4 shows cr xs as a function of the average flux / 
within the interval. A clear proportionality is seen with cr xs ~ 
(0.95 ± 0.10)/;. A log-normally shaped distribution with a pro- 
portionality of the excess variance with the flux are evidence of 
log-normal flux variations (Aitchison & Brown 1957). 

Below a flux level of few tenths of Crab, one might expect 
a low baseline flux level that can be described by a Gaussian. 
At these low fluxes, however, detector sensitivity threshold ef- 
fects become increasingly important, i.e. lower fluxes cannot 
be detected significantly within the short observation windows. 
Therefore, the mean of the Gaussian fit (0.33 Crab flux units) 
must be treated as an absolute upper limit on the integral base- 
line flux above 1 TeV from Mrk 42 1 . More information such as 




f / Crab 

Fig. 4. The excess variance <r xs as a function of the average flux within 
equal-length intervals. A line fit corresponding to <x„ ~ (0.95 ± 0. 10) ft 
is also given. 



the variability type (e.g. red noise / white noise / blue noise) or 
the duty-cycle of the object can be extracted from the flux state 
distributions. Previously, first steps in this direction were taken 
by Tluczykont et al. (2007); however, further investigations are 
necessary using simulations and will be the topic of a subsequent 
publication. 

With increasing exposure achieved by the major current gen- 
eration experiments more data will soon be available. As dis- 
cussed previously, one important aspect for long-term variabil- 
ity studies will be to carry out unbiased (random) observations 
to avoid systematic selection effects. 



4. Summary & outlook 

VHE y-ray flux measurements of AGN have been collected from 
observations of different experiments since 1992. For the first 
time these data were combined into single long-term lightcurves, 
and are provided for further analysis in the standard FITS file 
format. The BL Lac object Mrk 421 yields the most extended 
dataset with a combined lightcurve spanning more than a decade. 
The collected data is publicly available for download 2 

The observed flux states averaged over 12 h of Mrk 421 
and Mrk 501 above 1 TeV are consistent within statistical and 
systematic errors with a linear correlation with daily averaged 
RXTE/ ASM X-ray count rates, with correlation coefficients of 
0.65 and 0.68 respectively. 

The long-term flux state distribution of Mrk 421 can be well 
described by the sum of a Gaussian and a log-normal distribu- 
tion. This behavior can be interpreted by the observation of a 
baseline flux and a stochastic flux state population governed by 
an underlying multiplicative process. In the framework of this 
interpretation, the combined dataset on Mrk 421 allowed setting 
an absolute upper limit on the baseline flux of Mrk 421 of 0.33 
Crab flux units above 1 TeV. 

Unbiased data from monitoring campaigns by ongoing ex- 
periments (for a review see Weekes 2008) and data from future 
experiments such as CTA will further extend the long-term cov- 
erage of variable objects. A unified standard data format will 
become increasingly important in the future, when more high- 
sensitivity VHE data becomes available, and when interfacing 
with other fields of classical astronomy and astroparticle physics 
(multi-wavelength and multi-messenger) intensifies. We intend 



http : / /nuastro-zeuthen . desy . de/magicjexperiment/ 
pro jects/lightjcurve_archive/index_eng. html 
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to Store these future data as described in this work in Standard Weekes, T. C. 2008, in American Institute of Physics Conference Series, Vol. 
data formats using the introduced SLF field Structure definition. American Institute of Physics Conference Series, ed. F. A. Aharonian, 

W. Hofmann, & F. Rieger, 3-17 
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